Numerical stability of the Z4c formulation of general relativity 
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We study numerical stability of different approaches to the discretization of a conformal decom- 
position of the Z4 formulation of general relativity. We demonstrate that in the linear, constant 
coefficient regime a novel discretization for tensors is formally numerically stable with a method 
of lines time-integrator. We then perform a full set of "app/es with apples" tests on the non-linear 
system, and thus present numerical evidence that both the new and standard discretizations are, 
in some sense, numerically stable in the non-linear regime. The results of the Z4c numerical tests 
are compared with those of BSSNOK evolutions. We typically do not employ the Z4c constraint 
damping scheme and find that in the robust stability and gauge wave tests the Z4c evolutions re- 
sult in lower constraint violation at the same resolution as the BSSNOK evolutions. In the gauge 
wave tests we find that the Z4c evolutions maintain the desired convergence factor over many more 
light-crossing times than the BSSNOK tests. The difference in the remaining tests is marginal. 
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I. INTRODUCTION 

There are currently two main formulations used in the 
numerical evolution of astrophysically interesting space- 
times by the methods of numerical relativity. The first 
was pioneered by Pretorius [l|, and later adopted by 



other numerical relativity groups, for example [3|-|5|. In 
this approach the generalized harmonic gauge formula- 
tion of general relativity, GHG d, 01 1 is employed with 
black hole excision. That is, inside the black hole hori- 
zon the numerical mesh contains a "cut-out" region. The 
generalized harmonic formulation has a number of de- 
sirable properties. The first is that it has a trivially 
wave-like principal part, which allows the construction 
of boundary conditions that lead to a well-posed initial 
boundary value problem that can be implemented numer- 
ically [sl-lsl-fl^. That the wave-like nature of the formula- 
tion is inherited by the constraint subsystem means that 
these boundary conditions can conveniently be made con- 
straint preserving, and that when numerical error causes 
violations of the constraints, this violation may propa- 
gate off of the numerical grid, hopefully to be harmlessly 
absorbed by the aforementioned boundary conditions. 
The evolution system also admits a constraint damping 
scheme |15| . which has proven important in numerical 
applications to avoid constraint violating blow-ups. The 
second main formulation is the combination of the BSS- 
NOK system, the moving puncture family of g au g e condi- 
tions, and puncture black- hole initial data [1^113 ■ Using 
this method coordinate singularities, or punctures, are 
explicitly advected across the numerical mesh. The punc- 
ture method has proven extremely robust in the evolution 
of even extreme initial data [l^, [l^ . But on the other 
hand, whilst well-posedness for the initial value prob- 
lem of BSSNOK with suitable members of the puncture 
gauge family has been established [13] the complicated 
structure of the principal part has prevented the same 
development of boundary conditions for the system, al- 
though interesting progress has been made [2l| . 

This status begs the question: is it possible to write 
down a formulation with the strengths of both GHG and 
BSSNOK for numerical applications? In |22| a natural 
candidate, namely a conformal decomposition of the Z4 
formulation |23l - [29| . Z4c was identified, and indeed found 
to give favourable results over those of BSSNOK in the 
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context of spherical symmetry. Like that of GHG, the 
constraint subsystem of the Z4 formulation has a triv- 
ial wave-like principal part. This structure is inherited 
by the constraint subsystem of the conformally decom- 
posed system provided that the additional constraints 
introduced by the decomposition are explicitly imposed. 
This property allows the convenient construction of con- 
strai nt p reserving boundary conditions for the Z4c sys- 
tem [30[. The Z4c system furthermore inherits the con- 
straint damping scheme of Z4, which was studied in de- 
tail in [3l| . So far the numerical studies of Z4c have been 
restricted to spherical symmetry. Recently a variant of 
the conformal decomposition, CCZ4, was performed, and 
for the first time numerical evolutions in three dimen- 
sions were presented [32|. The difference between Z4c 
and CCZ4 is that in the conformal decomposition of Z4c 
non-principal constraint additions arc discarded in such 
a way to make the resulting evolution system as close 
as possible to BSSNOK. On the other hand in CCZ4 no 
constraint terms are discarded. Therefore the two formu- 
lations share the same principal part, and thus the same 
basic PDE properties. Now that there is robust evidence 
in favor of the Z4c system in spherical symmetry we also 
turn our attention to three spatial dimensions. Since 
there is a wealth of experience concerning the robustness 
of the GHG and BSSNOK formulations in applications, 
significant evidence must be presented that an alterna- 
tive formulation is competitive, and we therefore take a 
conservative approach, and focus here on the question of 
numerical stability. 

In section [n] we present the equations of motion of the 
Z4c system and the gauge conditions that we employ in 
both our analytical and numerical studies. In section Hill 
we introduce a novel discretization of tensors and show 
that in the linear constant coefficient approximation Z4c, 
coupled to the puncture gauge, is numerically stable with 
this discretization and a method of lines time-integrator. 
In section ITVl we present a complete set of apples with 
apples'^ tests for the Z4c formulation. In appendices |X] 
and [B] we discuss hypcrbolicity of the conformal decom- 
position without imposition of the algebraic constraints, 
and symmetric hyperbolicity of the formulation with the 
puncture gauge. We conclude in section fVl 

II. THE Z4C FORMULATION 

The Z4 formulation jl^, with constraint damping [T5| 
replaces the Einstein field equations Gab = 87rTaf, by 

Gab = 87rTab - V aZb — ^b^a + Qab^ cZ'^ 

+ Ki[naZb + UbZa + K29abncZ''], (1) 

where Gab ~ Rab~\Rgab is the Einstein tensor while Rab 
is the Ricci tensor, Va is the covariant derivative oper- 
ator compatible with gab and Za is an additional vector 
field of constraints. Under the standard 3-1-1 decom- 
position, against a timelike unit normal vector n° , defin- 
ing 8 = —n°^Za and Zi =-L° Za one finds the expressions 



given in [15|. If the spacetime is without boundary and 
does not admit a Killing vector, it can be shown that 
if the constraints are satisfied in one spacelike slice then 
they will all vanish all times H^l ■ We may therefore take 
a free-evolution approach to the problem. For PDEs and 
numerical analysis we work for in the expanded phase 
space, in which the constraints may be violated. In nu- 
merical applications the constraints are to be solved for 
initial data, and their compatibility with the evolution 
equations means that any violation should converge away 
with resolution. 

From the PDEs point of view, only the principal 
derivative additions can affect well-posedness of the ini- 
tial value problem of the formulation as a PDE system. 
In order to make a conformal decomposition we there- 
fore discard non-principal terms in such a way that the 
resulting equations of motion will have a form similar to 
those of the BSSNOK formulation when written in terms 
of conformal variables. Of course we keep the constraint 
dam ping terms. The Z4c formulation equations of mo- 
tion [22l.[30j are 

dtli] = -2aKij + Cpjij, (2) 
dfK^j = -DiDja + a[Rij + KK,^ - 2K,kK^ j + 2L)(,Zj) 
- Ki(l -f K2)07ij] 4- ATTa['^ij{S - Padm) - 25"^^] 
+ CpK,j, (3) 

for the metric and extrinsic curvature, and 

dtQ = ^a[F + 2b' Z, - 2ki(2 + k^)^] + CpQ, (4) 

+ P'b,z„ (5) 

for the constraints Q and Z^, where here we define 

b^Zj^r'^kAb~^'Z% (6) 

and use the shorthands 

H = R- K,jK'^ + _ I6^p^„„, (7) 
M, = [K,j - 7,, K] - 8^5,, (8) 

for the Hamiltonian and momentum constraints. Similar 
to the conformal transformation made in the BSSNOK 
formalism, we introduce conformal variables 

X = l~^'\ (9) 

= - \l^,K), k = K-2Q, (10) 

r = 27^^Z, + 7^^7"7jfe,/, f = f^fe7'^ (11) 
Under this change of variables the evolution equations 
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become 

dtX = 



-^XHK + 2Q)-D,(i\ 



(12) 
(13) 



-DiD'a + a[A,jA'^ + -{K + 26)^ 

O 



dt K ■ 



+ Ki (1 - K2)e] + 4^a[^ + Padm] + 13' m, (14) 



2A,fei%] + 2i,(,9,)/3'= 



(15) 



5*9 



2Ki(2 + K2)e] +p'd,e, 



(16) 



2a[r,feiJ"^ - ^A'^d, Inx - Jr^5,(2/^ + 9) 



(17) 



where we absorb the constraint addition into the Ricci 
tensor according to 



Ri 



(18) 
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4x 



i?,; 



(19) 



(20) 



In this study we consider only vacuum spacetimes, so 
Sij , Si and Padm all vanish. Note that the conformal 
decomposition introduces two algebraic constraints 



D = lndct7 = 0, 



T EE Y^A,, - 0, 



(21) 



to the system. In our numerical experiments they are im- 
posed explicitly after every time-step. If the constraint 
projection is not performed then one is performing free- 
evolution in the larger phase space in which the algebraic 
constraints are violated, meaning that the PDE proper- 
ties of the system must be re-evaluated. We close the 
system with the puncture gauge condition, which con- 
sists of the Bona-Masso lapse [S^ and the gamma-driver 
shift conditions ISJl, 



dta 



(22) 
(23) 



When coupled to the puncture gauge the Z4c formula- 
tion forms a PDE system that is generically strongly hy- 
perbolic. In our numerical applications we almost always 
take the "1+log" variant of the lapse condition = 2/a 
and /.is = 1/0?, although some tests are performed with 
harmonic lapse /iL = 1. The shift damping term 77 is 
chosen according to the needs of the numerical test. We 
also sometimes employ the harmonic shift condition 



dtp' = o? [xf' + \xTx., - Xl'^d, ln«] -I- p^d.P' 



(24) 



III. NUMERICAL STABILITY 



In this section we demonstrate that when coupled to 
the puncture gauge and linearized around the Minkowski 
spacetime, the Z4 formulation is numerically stable with 
a particular discretization. The calculations build on 
studies of numerical stability for first order in time, sec- 
ond order in space systems j35l - l37j . and essentially rely 
on the theory of (ssl. |39|. 



A. A novel discretization of second order systems 
of tensors 

Motivation: For linear, constant coefficient first or- 
der strongly hyperbolic systems, strong hyperbolicity is 
enough to guarantee numerical stability under a method 
of lines approach, which is no longer the case for first or- 
der in time, second order in space systems. If in addition 
to being strongly hyperbolic, the PDE system is symmet- 
ric hyperbolic, then a discrete energy method can be used 
to guarantee numerical stability. There are many second 
order PDE systems of interest that are only strongly hy- 
perbolic, for which other methods must be used in anal- 
ysis. Unfortunately when using the standard discretiza- 
tion in space for strongly hyperbolic second order systems 
it is often not even possible to analyze whether or not the 
resulting semi-discrete system is stable with computer al- 
gebra. The reason for this is that there is not a straight- 
forward relationship between the principal symbol of the 
semi-discrete system and that of the continuum; various 
sectors of the principal symbol that are decoupled at the 
continuum become entangled in the semi-discrete system. 
Here we present a novel discretization for systems of ten- 
sors for which semi-discrete stability analysis can at least 
be tackled. Intuitively this is made possible because af- 
ter Fourier-transform the novel discretization has only 
one wave- vector, whereas the standard discretization has 
two. The novel discretization does not guarantee numer- 
ical stability whenever the continuum system is strongly 
hyperbolic. Therefore the situation is still not entirely 
satisfactory. 

The standard discretization: For a grid- function /, 
where here and in the following we suppress indices label- 
ing the position on the grid, the standard second order 
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discretization accurate for second order in space systems 
is 



D, 



Oi, 



0^0, ^"^^"^ ^'^'^ (25) 



where we use the standard notation for centered and one- 
sided finite difference operators. The standard fourth 
order accurate discretization for second order in space 
systems is 



D 



(4) 



Do^ 1 



-D+,D^, , 



^(4)^(4) 



(26) 



if « 7^ 3, 



D+,D^,{1 ~ ^D+,D^,) if; 



3, 



where here we drop the summation convention. This 
approach to discretization can of course be naturahy ex- 
tended to higher order accuracy. 

A novel approach to the discretization of tensors: We 
are concerned with numerical approximations to the spa- 
tial derivatives of scalars, vectors and symmetric two- 
tensors, and assume that we have a metric r]ij in space. 
Suppose that we are given the finite difference approxi- 



mations and D\j' to the first and second derivatives, 
respectively. For the first derivative we must therefore al- 
ways use dI^\ For the second derivatives we distinguish 
as follows. The Laplace operator is approximated by 



l(2) 



ifd.dj A(2) EE ifDlf. (27) 
Second gradients of scalars arc approximated by 

d,d,u ^ [dI'^D^Y^ + l^j^^^^u. (28) 

For vectors we distinguish gradient-divergence terms by 
using repeated application of the first derivative, 



(29) 



Likewise for symmetric two-tensors we approximate 
divergence- divergence terms by repeated application of 
the first derivative approximation, although such terms 
do not appear in our applications. In [35| numerical 
stability of a parametrized generalization of the KWB 
formulation [40| of electromagnetism was discussed. It 
was found that with the standard discretization there are 
choices of the formulation parameter that render the sys- 
tem numerically unstable, despite the initial value prob- 
lem being well-posed at the continuum. It is straightfor- 
ward to show that with the novel discretization, using the 
standard second order discretization for the raw approx- 
imation, the semi-discrete system is numerically stable 
for all choices of the parameters that render the initial 
value problem well-posed. 



B. Application to the Z4 formulation 

The linearized system: In this section we present the 
Z4c field equations linearized around the line element 



ds^ = ~At^ + ^jdsMx-', 



(30) 



where we denote the constant background spatial metric 
by rjij. The background shift, extrinsic curvature and Z4 
constraints are taken to vanish. We denote the pertur- 
bations by a, 7ij, O, Zi^ Kij. The linearized equations 
of motion are 



9t7,, =-2K„- + 2%/3j), 



(31) 
(32) 
(33) 



dtKij = -^d'^dkjij - ^didj'-f - didja + d^Jjj, (34) 



(35) 



dtP = d'djP' + ^d,d,l3^ - p'K + 2d'Q, (36) 

where here we have defined jls = ('7)^MSj ^-nd used the 
Nagy-Ortiz-Reula (NOR) |4l| variable 

f = 2Z^ + r,^^f,''^[dkji,-ld,jki]. (37) 

The linearized Hamiltonian and momentum constraints 
are 



(38) 
(39) 



The semi-discrete system: We now discretize the sys- 
tem in space but leave time continuous. Using the novel 
discretization described in section fill Al we write the sys- 
tem as 

^a^J = -2/^,,+2i?(i)(,;/3,), (40) 
dta^-fiLiK~2&), (41) 
dtP' = tisP, (42) 

[[i?(i),:i?(i),ff+ir),,A(2)]« + i?(i)(J,), (43) 



(44) 



3 3 

(45) 

where we suppress indices labeling the grid, which we 
take to be 



Xt = (xii, 2/12,^13) = (11/1,12/1,^3/1), 



(46) 
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with ii = 0...iV— 1, h ^ ^ the spatial resolution, and 
restrict our attention to 27r-periodic solutions. Note that 
in contrast to the continuum system the constraint sub- 
system does not close for the semi-discrete system. It is 
not obvious how to achieve closure of the constraint sub- 
system at the semi-discrete level without using a "(£)o)^" 
type discretization. Unfortunately, as discussed in [35| 
the Dq norm is not strong enough for analytic consider- 
ations, and therefore do not consider it further here. 

Psuedo- discrete reduction to first order: Performing 
a discrete fourier transform in space, defining 



(47) 



and = -TT + 2ti/N, -tt + An/N, . . . , -|-7r, for i = 1, 2, 3, 
and then introducing the reduction variables 

e = ifl^a, = ifl^/3'\ tij = ifl^jij, (48) 

results in a large set of ODEs which can be written as 

2s 



dtjij 



-2K,. 



for the metric components. 



dt~e = -in+HL{k -2Q), 
dt4>' = in+fisf\ 



(49) 

(50) 
(51) 



(52) 
(53) 
(54) 



for the reduction variables, and finally 

1 



dtQ 
dtf 



■[s'[s 



Itf 



21 e- 



2n+''' 



1 



[s [siSj] + -riijn \l+iss^Jj), 



L + -isf 

3^2+ 2 



I s 

3lL 



\iss'k 



2iss'e. 



(55) 
(56) 
(57) 



for the remaining quantities. Here we have defined 



Si — —iD, 



(1) 



SiS , 



(58) 



and use an index s to denote contraction with the unit 
wave- vector s'. Whilst estimating the growth of the solu- 
tions we may discard the non-principal terms because for 
the standard second and fourth order accurate discretiza- 
tions the non-principal additions are bounded (Theorem 
5.1.2 [13 )• Note that for both the standard second and 



fourth order accurate discretizations we have the proper- 
ties. 



1 ^ 



(59) 



for some positive constants a, 6, which we will use in 
the following. Upto non-principal terms the evolution 
of the metric components a, "fij is trivial, so they may 
be discarded in the following discussion. Defining the 
unit wave- vector Si, and using it with fjij as a metric in 
Fourier-space to perform a 2 -I- 1 decomposition reveals 
that, under the novel discretization the principal symbol 
of the pseudo-discrete reduction decouples into a scalar 

{Lss,iqq,eJ',kss,kqq,Q,(t>''), vector {TsA, , K^A, 

and tensor (iASi-^Js) sector, just like in the continuum 
PDEs analysis. 

Characteristic variables: Each of the three sectors 
gcnerically admits a complete set of characteristic vari- 
ables. In the scalar sector they are 



—e ± ^/JUK, 



(60) 



which are associated with the lapse gauge condition, 
propagating with speed zL^/JIZi^l, and appear almost 
identically at the continuum level. For brevity we use 



the shorthands K 
Next we find 



Kss + K, 



11 



20 and A = t^s + ? 



qq- 



U±l 



2n 



3f72 



s2 s2 

— iss + 2n— 



± 12(r22 _ s^)Kss 



=F 6(0^ + 2s^)kgg ± 245^9 - 12ns/^ 



(61) 



propagating with speeds ±ifl, which corresponds to light 
speed in the continuum system. At the continuum four 
characteristic variables in the scalar sector propagate at 
light speed. In the semi-discrete system with our discrec- 
tization the third and fourth of these characteristics 



U±y_ 



3fis^^ — + s'^{p,s - 1) 
17+(w2 _^^r!2) 

4s2(/is - 1) 



A 



±2 



(ml - Ms)(3w_ 



4n^)^s^^iL{fls-l) 



fisK 



02 



Ms - 1 



(62) 



become coupled to the characteristic variables of the 
shift. They propagate with speeds ±u_, where we de- 
fine v± by 



4 



1 



± l\J [(4 + SfisW + (MS - - 48/2sf^4. (63) 



Under the substitution J7 — s the continuum result (in 
fourier space) is recovered. The last of the characteris- 
tic variables in the scalar sector are naturally associated 
with the gamma driver condition. They are given by the 
same expression as (|62p under the replacement v- — >■ , 
and have speeds ±w-(. . Special cases in the choice of the 
gauge parameters that prevent the existence of a com- 
plete set of characteristic variables in the scalar sector 
are discussed below. Thankfully in the vector sector we 
have the straightforward 



n 



(64) 
(65) 



with speeds ±zf2 and Finally in the discretiza- 

tion of the gravitational wave degrees of freedom we have 
the trivial, 



TF 

'■AB ±1 



ml 



''AB =■= ^ t 



^AB^ 



(66) 



with speeds ±i51, and where TF denotes that the trace is 
removed with respect to the projected background met- 
ric. 

Special cases: Assuming that /^/, > 1, there are two 
special cases in the choice of the gauge parameters fiL , fis 
under which the scalar sector does not admit the com- 
plete set of characteristic variables (|60ll62p given above. 
They are 



[(4 + 3fls)n^ - 
(3/iL-4)(AiL 



(Ms-l)s']'-48Msf^'* = 0, (67) 
- /is)f^'- ML (MS -l)s' = 0, (68) 



which correspond to Ap^s ~ 3q!^ and 3/1^ = 4/2^ at the 
continuum in curved space. These special cases are more 
subtle in the semi-discrete system because instability de- 
pends on both the background metric and the specific 
order of the discretization. Therefore we restrict our dis- 
cussion here to rjij = Sij, the gauge parameter ml = 2. 
Then for either /2s = | or /X5 = |, in the limit of high 
resolution, /i — >■ the characteristic variables become de- 
generate in a neighborhood of the grid mode = 0. 
This situation is the same as that of the ADM formula- 
tion discussed in [s^ . Any other mode satisfying (|67II68|) 
will be unstable. There are no such modes for /is < | 
or I < Ms < I because for the standard second and 
fourth order discretizations a = 1 in (|59p above with the 
Minkowski background. In the limit of high resolution 
one also expects to find unstable modes in the region 
I < Ms < 2, although for particular choices of Ms they 
may also appear at low resolutions. The general case is 
sketched in Fig. [1] 

Numerical Stability: Notice that at the lowest and 
highest grid-frequencies the principal symbol takes a dif- 
ferent form, because cither or vanish. As discussed 
in [3^ this poses no difficulty in the constant coefficient 




FIG. 1: The unstable gridmodes with the novel discretization 
and hl ~ 2. The left (blue) line denotes (|67|) and the right 
(red) (|68p . The dashed vertical lines depict the choices of Ms 
that render the system weakly hyperbolic at the continuum 
level. For those choices no stable discretization is possible. In 
practice for the weakly hyperbolic choices of p.s the character- 
istic variables become singular at high resolution, preventing 
the construction of a symmetrizer independent of h. In princi- 
ple one expects to be able to find a stable discretization for all 
other values of fls, but our novel discretization does not fulfil 
that requirement. The unstable modes should be truncated in 
the y-axis at some finite a > 1 determined by the discretiza- 
tion and background metric, and so large values of unsta- 
ble s/Q are not to be considered troublesome. On the other 
hand behavior of the red line in the region 3/iz,/4 < Ms < Mi 
is undesirable, because we always expect to find modes in 
that region on the numerical grid. Note that for our standard 
choice, fls = 1, the blue line diverges, provided that 77 ~ 1, 
implying that our standard choice is stable, at least in that ap- 
proximation. There are never unstable modes when iil < Ms • 



case, provided that the principal symbol remains diago- 
nalizablc for those frequencies, as the full solution may 
be written as a direct sum over the modes. The prin- 
cipal symbol vanishes at the lowest frequency and so is 
trivially diagonalizable. At the highest frequency, by con- 
struction, the principal symbol can be written as that of 
a system of decoupled wave equations, and so is also triv- 
ially diagonalizable. The argument given in (35| guaran- 
tees that standard method of lines time-integrators will 
not violate the reduction constraints. The existence of 
a complete set of characteristic variables for the system 
for every grid frequency guarantees the existence of a 
symmetrizer H{si), which in turn guarantees numerical 
stability, namely that the estimate 



Mtn,-)\\hM+ </^e"*"||u(0,-)||/v 



(69) 



holds for sufficiently small h, where a,K are constants, 
independent of h, that are not to be confused with the 
lapse and extrinsic curvature, provided that the Courant 
factor A, the ratio of the time-step and the spatial resolu- 
tion, is chosen sufficiently small. The Courant condition 
is A < q;o/(2v^3ml), in the flat-space, second order accu- 
rate case with Ms = 1 a-nd the stability radius of the 
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time integrator. The norm || • ||h,_D^ is given by 

3 3 



2—1 i-.j — ^ 

3 3 3 

+ Eii/'ii' + E iii^.,ii^+Eii^+HP. 



3 3 

E \\D+^|3'\\l+ E 

— l i,j,k—l 



(70) 



Discussion: The inclusion of a constant background 
shift is straightforward since it affects the characteristic 
speeds in a trivial way, and has no affect on the character- 
istic variables. Different choices of discretization of the 
shift advection derivatives were studied in detail in [36j . 
Setting the background lapse to unity was done only for 
convenience and does not affect the results qualitatively. 
In the linear constant coefficient regime the conformal 
variables can be expressed as a linear combination of the 
NOR variables that we have used in our calculations. The 
transformation to conformal variables is not however just 
a change of variables, because the transformation intro- 
duces the algebraic constraints D and T. In the linear 



constant coefficient discussed in [37[, if an ex- 

plicit polynomial time-integrator is used then linearity 
of the system guarantees that the algebraic constraints 
will not be violated if they are initially satisfied. The 
time-integrator can furthermore be modified by a con- 
straint projection step, so that even if the initial data 
does include violate the algebraic constraints then nu- 
merical stability of the algorithm is trivially recovered. 
In the variable coefficient and non-linear cases, the con- 
straint projection step is essential because error will oth- 
erwise violate the algebraic constraints, and as shown in 
Appendix [XI when the algebraic constraints are violated 
the Z4c formulation is generically only weakly hyperbolic. 
It would be desirable to extend the results presented in 
this section to the variable coefficient case, rather than 
linearizing and freezing coefficients. The most power- 
ful method for doing so would be to use a summation 
by parts finite-differencing scheme to conserve a semi- 
discrete energy. Unfortunately such an approach does 
not work because, as shown in Appendix |Bl even at the 
continuum level no such energy exists. 



IV. THE APPLES WITH APPLES TESTS 

A. Numerical setup 

For our numerical tests, we use the AMSS-NCKU 
code [i^ - li^ . The code employs MPI parallelization for 
moving box mesh refinement and method of lines time in- 
tegration. In physics applications the code uses the mov- 
ing puncture gauge with BSSNOK to evolve black hole 



spacetimes. Here it is used in a much simpler context: 
each of the Apples with Apples tests employs just a single 
mesh, and peridodic boundary conditions in space, which 
avoids implementing complicated conditions in 3D. Fol- 
lowing [45I, we set simulation domain x £ [—0.5,0.5]. 
In contrast to the normal AMSS-NCKU setup, we use an 
unstaggered, or vertex centered, grid a;„ = —0.5 + ndx, 
n = 0...50p. Wc take dx — dy ~ dz ~ l/(50p), 
p ~ 1, 2, 4, and so on. In the y and z directions, fol- 
lowing [4^, we use only five grid points. The Courant 
factor A = dt/dx is taken to be either 0.1, 0.25, or 0.5 
according to the specific goal of each test. For spatial 
derivatives, we restrict exclusively to second order accu- 
rate finite differencing. For advection terms /3'(?i, we use 
the upwind scheme 



^(up) ^ I D 



ldx,D+,D+,ii 
+ ldx,D^,D^, li < 0. ^ > 



For the remaining spatial derivatives we employ either 
the standard spatial discretization (j25j) . or the novel 
scheme described in section IIII Al In |47| it was shown, 
albeit for gauge conditions other than those that we em- 
ploy here, that without artificial dissipation the BSSNOK 
formulation is not numerically stable. Since we are com- 
paring the Z4c formulation with BSSNOK, we therefore 
always use Kreiss-Oliger numerical dissipation [s^ 



dtu~^dtu-(j^dxl{D,+ D,-)' 



(72) 



with dissipation parameter cr = 0.02. In the original Ap- 
ples with Apples tests, gauge conditions were prescribed. 
Unfortunately some of those prescriptions result in an ill- 
posed initial value problem when coupled to either the 
BSSNOK or Z4c formulations. Furthermore we prefer to 
study error properties of the exact system that we will 
later use to perform astrophysical simulations. We there- 
fore use only the puncture gauge (|22II23|) with fi^ = 2/a, 
fis = ^jo? and rj = 2 unless stated otherwise. We are 
mainly concerned with the effect of each formulation on 
the accuracy of the numerical calculation at finite reso- 
lution and the effect of the novel discretization scheme 
versus the standard discretization for the Z4c system. In 
our tests we frequently use the constraint monitor 



" 1 ,/WT':i~~MmrT'^i~&& 



for Z4c, 
for BSSNOK, 
(73) 



to assses the quality of the numerical solution. Here H 
is the Hamiltonian constraint violation, M* is the mo- 
mentum constraint violation and G" is the Gamma con- 
straint violation. Since 2Z* corresponds to the Gamma 
constraint in BSSNOK formalism, we use A'jijZ'^Z^ to 
mimic the Gamma constraint violation part. 
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FIG. 2: Comparison of the constraint violation monitor (|73|) for the robust stability test with 77 = 2. The left subplot compares 
BSSNOK and Z4c, where the standard discretization scheme is used. Numerical experiments reveal that the drop in the 
constraint monitor in the Z4c tests is caused neither by artificial dissipation nor the t] parameter. The right subplot compares 
the standard discretization scheme and the novel discretization scheme for Z4c formalism. 



B. Robust stability test 

First we perform the robust stability test. Besides 
the gauge choice we follow the original presciption. Ini- 
tially all variables including lapse and shift are set to 
the Minkowski value perturbed by a random number 
e e {-lQ-^^/p^,lQ-^^/p^). We set the Courant factor 
to A = 0.5. In Fig. [51 the constraint violation monitor 
respect to time is plotted for evolutions with 77 = 2. In 
the left subplot we compare Z4c formalism with BSSNOK 
formalism using in each case the standard discretization 
scheme. Although neither BSSNOK nor Z4c suffer from 
rapid undesirable growth of the solution, the constraint 
violations arc slowly increasing for BSSNOK formalism 
while they are decreasing for the Z4c test. By the end 
of the evolution, the constraint violation in the Z4c tests 
are several orders of magnitude lower than those of the 
BSSNOK test. We do not have an explanation for this 
behavior. The usual argument is that the propagating 
constraint subsystem of the Z4c formulation should pre- 
vent constraint violation from growing on the grid, but 
in this test the constraint violation can not be absorbed 
by the boundary conditions. To investigate the cause 
of the decay of the constraints we have performed tests 
without artificial dissipation. There we find that the con- 
straint violation levels off around an order of magnitude 
above the value obtained with artificial dissipation. We 
also tried evolutions with r/ — 0, and find that the con- 
straint monitor is larger than with r/ = 2, but that the 
qualitative behavior is unaltered. Experimenting with 
the constraint damping scheme with ki = 0.02, K2 = 0, 
we find a minor improvement in the constraint monitor, 
in line with our expectation from the results of |3l| . In 
the right subplot, we compare the standard discretization 
scheme and the novel scheme, for which we have a proof 
of numerical stability in the constant coefficient approx- 




FIG. 3: Comparison snapshots of 7^^ — 1 at t = 1000 for Z4c 
formalism with the standard discretization scheme. 



imation, each with the Z4c formalism. Here however we 
would like to reiterate the important point made in [s^, 
namely that formal numerical stability neither implies, 
nor is implied by robust stability. The robust stability 
test, in the form studied here, should be viewed only as 
a gauge of the errors present in the numerical evolution 
at a finite resolution. In appendix |X] wc present more 
demanding convergence tests. 



C. Linear wave test 



The second test is the evolution of a linearized gravita- 
tional wave on the Minkowski background. We set initial 
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data as 

Ixx = 1, Jyy = 1 + &, Jzz = l~b, 



— dth, Kzz = 

2 2 

and the remaining variables to zero. Here we choose 



a = Kyy^-dA K,, = ~-dtb, (74) 



b ^ A sin 



2Tr{x - t) 



(75) 



and d — 1, A — 10^^. The Courant factor is cho- 
sen A = 0.5. Such a small perturbation of the metric 
ensures that the numerical evolution remains essentially 
in the linear regime. Here the initial data are constraint 
violating, and the puncture gauge condition is not nec- 
essarily compatible with simple advection of the wave 
profile. In practice however it appears that to a good 
approximation the solution is a simple travelling wave; 
in Fig. [3] we compare the final jyy numerical profile, ob- 
tained with Z4c, to the analytic solution of the linearized 
Einstein equations with unit lapse and zero shift. There 
it appears that at least the phase of the numerical solu- 
tion is converging with resolution to the analytic solution. 
A more meaningful evaluation of the errors in the system 
is obtained by performing self-convergence tests on the 
data, for which we use the norm ([70)1 evaluated on 
the conformal variables of the Z4c and BSSNOK formu- 
lations. There we find perfect second order convergence 
in the norms for both BSSNOK and Z4c (Fig. H BSS- 
NOK gives the same behavior as Z4c, not shown), and 
estimate that after 1000 light crossing times, the error is 
approximately 2 X 10"^ for resolution p ~ 4 and 8 (Fig. [4]). 
Since the initial data for these evolutions is constraint vi- 
olating we do not present the constraint monitor for this 
test, although we find that here it stays around 1 x 10~^^ 
with either formulation. The novel discretization for Z4c 
gives almost identical results to those of the standard 
discretization. We tested the effect of constraint damp- 
ing terms and different values of rj in the gamma driver 



shift. Neither has an effect on the result of these tests. 
Concerning coordinate gauge conditions, for comparison, 
we have performed additional evolutions with harmonic 
lapse, = 1 in with harmonic shift ([M)) . We find 
both choices give very similar same convergence behavior 
and errors. 

In addition to above sine wave test, we also ado pted 
the Gaussian shaped linear wave test suggested in [48| . 
Instead of Eq. ([TS)) , the b takes 



b = ^exp[- 



2w2 



(76) 



with A = 10"^. In order to get periodic profile initially 
we set w = 0.5 which makes b roughly at a; = ±0.5. 
Since the effective wavelength for this gaussian wave is 1 
which is one tenth of the wave length of above sine wave, 
we need larger resolution to get reasonably convergent 
result. Reducing numerical dispersion requires more res- 
olution still. In all we find resolutions p = 16, 32 and 64 
can produce good convergent result (Fig. [5]). The error 
for p — 6A and 32 is comparable to the p = 8 and 4 of 
sine wave case, which is roughly 2 x 10^^. 
As to the two dimensional sine wave test 



7x 



1 + — 



H 
'2 



A' 



H = Asm 



2tt{x -y- V2t) 



1 - H, 
2 



1, 



(77) 



d ^ 1, A ~ 10^^ and other variables zero, the result 
is similar to the one dimensional test. But because the 
effective wave length for two dimensional sine wave is 
l/-\/2, a little higher resolution (p = 2,4,8) is needed to 
get comparable convergence behavior as in the one di- 
mensional case (Fig. [S]) . As in the earlier robust stability 
tests the constraint damping scheme with ki = 0.02, K2 ~ 
reduces the constraint violation very little. 
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FIG. 5: Convergence factor of gauss shaped linear wave test with amplitude A = 10 Puncture gauge is used. 
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FIG. 6: Convergence factors and and errors for as in Fig. [5] but for the diagonal sine wave test. 



D. Gauge wave, and shifted gauge wave tests 

The third and fourth tests are the unshifted gauge wave 
and shifted gauge wave tests. For the gauge wave test we 
take 



K 



lyy = 1. 



Iz 



1, 



a = Vl -h, 



as initial data. For the shifted gauge wave we use 



(78) 



1XX = 1+ b, Jyy = 1, 'Jzz = 1, (79) 

dtb 1 b 



K. 



2VTTb' 



VTTb 



px 



l + b' 



and in either case set the the remaining variables to 
zero. We again choose (|75p . but now with d = 1. 
Both A = 0.01 and A = 0.1 are tested. We set the 
Courant factor A = 0.5 in the unshifted case and A = 0.25 



in the evolution with non-trivial initial shift. Otherwise 
the grid setup is the same as in the robust stability test. 
It is well-known that with Harmonic lapse, ^.^ — 1^ and 
vanishing shift the BSSNOK formulation suffers from 
large undesirable growth of the errors in both of these 
tests H^. We recover that behavior for BSSNOK, and 
find that with the Z4c system the blow-up is delayed. For 
a direct comparison with the gauge wave results of [32[ we 
performed this test with Harmonic lapse, vanishing shift 
and constraint damping and ki = 1,K2 =0. There, in 
agreement with (32j we find that the errors do not result 
in a code crash, although at this resolution the experi- 
mental convergence factor is very poor. In contrast to 
the Harmonic lapse gauge condition, the puncture gauge 
is able to evolve the gauge wave for 1000 light crossing 
times with no sign of blow-up without constraint damp- 
ing. Since both A = 0.01 and A = 0.1 give roughly the 
same result, in Fig. [7] we plot the convergence factor and 
the norm of the error only for A = 0.01. As the 
first test of nonlinear effects in "apples with apples" test 
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suite, Z4c gives better convergence behavior than BSS- 
NOK. As shown in the left subplot of Fig. [71 Z4c can ap- 
proximately maintains a convergence factor around two 
while BSSNOK drops to negative values. Moreover, Z4c 
gives more than one order smaller error as shown in the 
right subplot of Fig. [7l In terms of constraint violation 
we find that without Z4c constraint damping the con- 
straint monitor is comparable for the two formulations 
throughout the evolutions. On the other hand we find 
that with Ki = 0.02, K2 = 0.0, the Z4c constraint vi- 
olation is an order of magnitude smaller than that of 
BSSNOK at the end of the evolution. The experimental 
convergence factor is unaffected by the use of the con- 
straint damping scheme in this test. Furthermore we 
find that with the puncture gauge condition the coor- 
dinates are rapidly driven to match incrtial coordinates 
of Minkowski space. In others words, the amplitude of 
the gauge wave decreases, as shown in Fig. [H so in some 
sense, the puncture gauge condition possesses a symme- 
try seeking property. As shown in this figure, we find 
that other gauge conditions do not have this property and 
the nonzero value of rj is also essential for the successful 
evolution. Note that this gauge condition is similar to 
that used in binary black hole simulations. The novel 
discretization gives near identical results to the standard 
discretization. Tests with the harmonic gauge result in 
both larger absolute error and larger constraint viola- 
tions. 

As for the two dimensional linear wave test, we use 



"fxx — Jyy — 1 ^5 Ixy — ~^ , 



H 
dtH 



Izz = 1, 



-Kx^ = — Kzz = — dfH, 



(80) 



H = A sin 



2tt{x - y - V2t) 



a = Vl - H, 



with d = 1, A = 0.01,0.1 for initial data for the two 
dimensional gauge wave test, which corresponds to the 
wave propagating along the diagonal direction of x — y 
plane. The results are similar to the one dimensional 
tests. 

As for the shifted gauge wave test, original "apples 
with apples" test suggested A = 0.5. We find such a large 
amplitude results in code crashes for both formalisms 
with any coordinate gauge. Our findings are similar to 
those found with a pseudo-spectral method as used in the 
SpEC code [4^ . The unshifted gauge wave tests also fail 
with amplitude A = 0.5. Restricting to A = 0.1, 0.01, we 
find the resulting error and the convergence behavior is 
very similar to unshifted case. As an example we compare 
the two results with A = 0.01 in Fig. [9] 



E. Gowdy wave test 



The final tests are the expanding and collapsing Gowdy 
waves. Following the suggestion of [45|,|4^, we choose for 
initial data 



K^x = ±-t-^''e^l\t-^ - A,,), a = r^'^e'^'^ 



Kyy = ± V/4e-^/4e^(-l - tP^t), 



(81) 



where sign is for the expanding Gowdy wave test and 
'-' sign is for the collapsing Gowdy wave test, and other 
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FIG. 8: Gauge wave test with the puncture gauge condition (|22II23I 77 = 2 or 0), with 1+log slice and zero shift and with 
harmonic lapse and harmonic shift. Here 7j,j, — 1 is plotted at times 0, 2, 6 and 14 respectively, for a test performed with the 
Z4c formalism and the standard discretization scheme. With puncture gauge and r] = 2 the gauge wave is rapidly suppressed 
by the gauge condition and at later times 7yy — 1 tends to zero. With ?7 = the code crashes soon after t — 7.2. While with 
other gauge conditions we can not see this symmetry seeking property. 



variables zero. Here P = Jo(27ri) cos(27ra;) and 
A = -27rtJo(27rt)Ji(27rt)cos^(27rx) 

+ 2TTH^[J^{2TTt) + Ji (27rt)] 

- ^[(27r)'[Jo'(27r) + J^{2Tr)] - 27rJo(2^) Ji(2^)] (82) 

where J„ are Bessel functions. Initially we set t = 
9.8753205829098 following [ii,!!^ exactly. The Courant 
factor is set as C = 0.05. Here the use of the punc- 
ture gauge is troublesome because it does not reproduce 
the preferred analytical gauge in the evolution. Unfor- 
tunately simply choosing the lapse and shift apriori re- 
sults in an ill-posed initial value problem, highlighting 
the fact that in numerical relativity there are still open 
questions related to the choice of gauge conditions. Ob- 
viously conditions must be chosen that render the initial 
value problem well-posed, but on the other hand such 
conditions may not correspond to those best-suited to 
study a problem analytically. 



The test results of expanding Gowdy wave are plot- 
ted in Fig. [TUl Although this test is very tough (46j . 
upto 1000, neither the BSSNOK formalism nor the Z4c 
formalism suffer from large growth of the errors. Our 
results demonstrate the robustness of puncture gauge 
condition. In the left subplot, we show the behavior of 
relative error for two different resolutions where p = 2 
and p = 4 are used. From this result we can see the 
increasing behavior of the error roughly follows power 
function of t instead of the exponentially increasing re- 
sult obtained by the spectral code in [4^. In the right 
subplot we show the constraint monitor with respect to 
time. In Fig. [11] we compare the convergence behavior 
with respect to Hamiltonian constraint violation. We 
find that Z4c formalism achieves better convergence re- 
sults than BSSNOK formalism at this resolution. Once 
again the novel finite difference scheme gives almost ex- 
actly the same result as the standard discretization. If 
the damping term with ki ~ 0.02 and K2 = is used. 
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FIG. 9: Comparison of the D+ norm of the error for itp=2 — 
Up^i between unshifted gauge wave test and shifted gauge 



wave test. 



the constraint violation is significantly reduced, but, as 
shown in Fig. [TOl does not reduce the norm of the dif- 
ference between two resolutions. Unfortunately in Fig llll 
we see that the damping term does reduce the experi- 
mental convergence factor. 

For the collapsing Gowdy wave, in contrast to the pro- 
posal of original ^'apples with appled^ tests, we do not 
transform the time coordinate to avoid singularity by 
hand. Instead we simply change the sign of initial extrin- 
sic curvature Kij in Eq. (|8T|) and evolve with the puncture 
gauge. Until 1000 light crossing times, when we stop the 
test, the code does not crash. As the solution approaches 
the singularity, the puncture gauge makes the evolution 
effectively stop by the collapse of the lapse. Both BSS- 
NOK and Z4c formalism gives almost identical results. In 
Fig. [T2|we show several snapshots for the lapse and shift 
vector, and the singularity avoiding property of punc- 
ture gauge condition can be clearly seen. Since BSS- 
NOK and Z4c formalisms give the similar result, only 
BSSNOK results are shown. For the collapsing Gowdy 
test we find that the effect of the constraint damping 
scheme, with ki = 0.02, K2 = 0, is negligible. 



V. CONCLUSIONS 

In order to demonstrate that a conformal decomposi- 
tion of the Z4 formulation of general relativity may be 
a useful tool for numerical relativity simulations in three 
spatial dimensions, we have studied numerical stability of 
the system coupled to the moving puncture gauge condi- 
tion. First, we proved that in combination with a novel 
discretization scheme, and periodic boundary conditions 
in space, the Z4c system is numerically stable when lin- 
earized around a constant background metric. We would 
like to extend our results to the variable coefficient case, 
but since Z4 coupled to the puncture gauge appears to 
be only strongly, but not symmetric hyperbolic, see ap- 
pendix |B1 there are few methods to hand. We then per- 



formed a complete set of the "apples with apples^^ tests, 
which we modified slightly to take into account the sta- 
tus of the initial value problem from the PDEs point of 
view, and discussed for each the relative behavior of Z4c 
in comparison with BSSNOK. We summarize the results 
of the tests as follows: 

Robust stability: In the robust stability tests we found 
that the Z4c evolutions result in a smaller constraint vi- 
olation, a feature that was insensitive to the choice of 
gauge. 

Linear waves: In the linear wave test we found only 
marginal differences between the two formulations, re- 
gardless of gauge conditions. It is perhaps expected that 
in this test the differences will be small, because the 
propagation of gravitational waves is of course a feature 
shared by both systems. 

Gauge waves: For the gauge waves test we found, 
using the puncture gauge, that the Z4c evolutions ex- 
hibit an experimentally computed self-convergence factor 
closer to the expected second order than the BSSNOK 
data over many light crossing times. In this context we 
found that the puncture gauge has a strong symmetry 
seeking effect and rapidly pushes the lapse to a constant, 
and the shift to zero. 

Gowdy waves: To evolve expanding and collapsing 
Gowdy spacctimes we again employed the puncture 
gauge. Here we find that although neither Z4c or BSS- 
NOK results in code crashes, long-term convergence is 
very difficult to achieve. Another problematic feature 
here is that although the puncture gauge results in a 
well-posed initial value problems for each of the formula- 
tions, it does not track the preferred coordinates on the 
spacctime under consideration. 

In appendix we also studied the effect of the al- 
gebraic constraint projection on numerical convergence. 
When the algebraic constraints are violated, the Z4c for- 
mulation is only weakly hyperbolic, which is reflected in 
numerical approximation as a failure to converge to the 
continuum solution, even for trivial initial data such as 
noise on top of the Minkowski spacctime. 

The evidence for the usefulness of a conformal decom- 
position of the Z4 formulation in astrophysical applica- 
tions is mounting. In spherical symmetry, the system was 
shown to have favourable properties over BSSNOK [2^ . 
especially in the evolution of spacetimes containing mat- 
ter. The trivial wave-like nature of the constraint subsys- 
tem was then used to construct high-order constraint pre- 
serving boundary conditions [lOl, which were studied in 
numerical applications in spherical symmetry. The con- 
straint damping scheme for the formulation was studied 
in detail in |3l| . There, once again in spherical symme- 
try, the efficacy of the constraint damping scheme was 
studied in numerical applications and it was found that 
the constraint damping scheme is a useful tool for the 
suppression of violations in vacuum spacetimes, but that 
when matter is present more work is required to paint a 
clear picture. For the first time, numerical evolutions of 
binary black-hole spacetimes with a conformal dccompo- 
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FIG. 10: Expanding Gowdy wave test with our puncture gauge condition. Courant factor C = 0.05 is used. Left subplot shows 
the D+ norm of the sohition difference between resohition p = 4 and p = 2. Right subplot compares the constraint monitor 
respect to time for BSSNOK formalism and Z4c formalism. For Z4c formalism we only show novel finite difference scheme here 
because standard finite difference scheme gives almost the same result. Here for Z4c formalism, with (ki = 0.02, K2 = 0) and 
without damping term are also compared. 
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FIG. 11: Expanding Gowdy wave test with our puncture 
gauge condition. Courant factor C — 0.05 is used. Conver- 
gence factor log2 ( ^^^^2) )• Here H{p) means the Hamiltonian 
constraint violation for resolution p. 



boundary value problem and can be implemented in a 
3D production code. Second, conclusive evidence of the 
benefits of Z4c over BSSNOK is needed in applications. 
We hope to address both of these issues shortly. 
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Appendix A: Z4c without enforcement of the 
algebraic constraints 



sition of Z4, CCZ4, were presented in [32|. Here we have 
focused on formal numerical stability, and compared nu- 
merical results obtained with BSSNOK and Z4c in a sim- 
ple context. Overall we find some benefit to the use of 
Z4c, although in some sense the tests we have performed 
here, namely the evolution of vacuum spacetimes with 
periodic boundary conditions, arc not optimally suited 
to the advantages found in spherical symmetry. In our 
view two points remain to be addressed before one can 
consider a whole-sale replacement of BSSNOK in applica- 
tions. First, it is desirable to obtain boundary conditions 
that are constraint preserving, minimise the incoming 
gravitational wave content, lead to a wcU-poscd initial 



1. Weak hyperbolicity with the puncture gauge 

Without assuming that the algebraic constraints are 
enforced, the equations of motion for the Z4c formula- 
tion can not naturally be written in terms of the ADM 
variables. Under the standard 2 -t- 1 decomposition [i^ 
and upto transverse and non-principal derivatives, they 
are given by 



dtX 
dass ■ 



;XHK + 28) - + P'dsX, 



-2aAs 
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FIG. 12: Collapsing Gowdy wave test. Courant factor is 0.05. BSSNOK formalism is used with resolution p — 1. Left subplot 
is the snapshot of a at different time; right subplot is for a approaches to freeze the evolution. increases at very 
early time to adapt the evolution and decrease to also soon later. They together make the evolution freeze near the universe 
singularity. 



dt%q = -2aAgg - -xdsP' + P'ds%q, (A3) 



dta = —fi^a^K + P^dgOt, 



(A4) 
(A5) 
(A6) 



dtAss = —dsdsX ~ —dsds^ss + —dsds%q - -xdsdga 



(A7) 



~ Ct cx 2 

dtAqq = --dsdsX + i^dsds^ss - -^dsdslqq + ^xdsdsOi 
6 6 6 3 



-ax^dst' +P'dsAqq, 



dt& = -dsdsX - —dsds^ss - —dsds^qq + ^rX^^r 
X 4x 4x 2 

dtT' = —dsdsp' - —d.sK - —9,9 + P'dsT' 
3x 3x 3x 



(A8) 
(A9) 



(AlO) 



in the tensor sector. Here the analysis is essentially for 
the linearized system, and we use the background met- 
ric to raise and lower indices. We also use the con- 
vention that the spatial vector s' is normalized to one 
against the background physical metric rather than us- 
ing the conformal metric. The principal symbol of each 
sector can be trivially read off from these equations. In 
geometric units the tensor sector has speeds ±1, and a 
complete set of characteristic variables. Likewise the vec- 
tor sector has speeds ±(1,^^^), with fls = Ms/x and a 
complete set of characteristic variables. The scalar sec- 
tor has speeds 0, ±(1, y/JIE, yAflsJ^), but does not have 
a full set of characteristic variables; the formulation is 
thus only weakly hyperbolic, and so does not have a well- 
posed initial value problem. It may be possible to modify 
the equations of motion by adding (possibly derivatives) 
of the algebraic constraints D and T to the equations of 
motion to achieve strong, or even symmetric hyperbolic- 
ity. But in any case it is not true that when the algebraic 
constraints are violated the PDE properties of the system 
arc unaltered. 



in the scalar sector, 

dtlsA = -2aAsA + P^'ds^sA, 



(All) 
(A12) 

dtA.A = -^dsdsisA + ^xc'.f ^ -f P'd,A,A. (A13) 
2x 2 



X 

in the vector sector and 



(A14) 



(A15) 



daX% = -2aAX%+[3^ds^XB, 
dtAYs = -f 9,a,7ll + P'dsIZ, (A16) 



2. Numerical convergence without constraint 
projection 

In the previous subsection we saw that when the alge- 
braic constraints arc not enforced, the Z4c formulation 
is only weakly hyperbolic. The open question is how ill- 
posedness of the initial value problem will be reflected 
in numerical approximation. To answer the question we 
follow [131 and perform convergence tests, using the D+ 
norm, as defined by (|70p . using initial data inspired by 
the robust stability test. Specifically we set the initial 
perturbation to e e (— 10~"^/p^, 10~^/pP). For the metric 
components we set p — 3, and use p = 2 for the remain- 
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FIG. 13: Convergence test with the puncture gauge condition (|22)I23|) . Convergence factors (|A17|) in time for different resolution 
are shown. The left subplots are without algebra constraint enforcement while the right ones with. Here the adjusted finite 
difference scheme is used. 



ing variables. This choice is needed to guarantee sec- 
ond order consistency of the initial data in the D+ norm, 
which is of course necessary for convergence in that norm. 
We employ the same numerical grid used in the robust 
stability tests, but instead of evolving for of 1000 light 
crossing times, we consider the interval t £ [0,0.8]. We 
set the Courant factor to A = 0.5, and use artificial dis- 
sipation ([72]) as usual with a = 0.02. We consider higher 
resolutions than those used in the robust stabiliy test, 
and compute additionally data with p = 16, . . . , 1024, in 
factors of two. Before calculating the the Z?+ norm, we 
restrict the solution to the p = 1 grid. The convergence 
factor with respect to resolutions p and 2p is 



1 f I K^P '^^cxact 1 1 
= loan r—^ 



C = log 



\\\u2p 



Ucxact Da 



(A17) 



The convergence factors obtained with the puncture 
gauge condition (|22j|23p and the novel discretization are 
shown in Fig. 1131 Very similar results are obtained with 
the standard discrctiazation and we interpret the close 
agreement of the two tests as numerical evidence that 
the standard discretization is formally numerically stable, 
with constraint projection, despite our inablility to prove 
so. At the lowest three resolutions the constraint projec- 
tion does not seem to have any effect on the result, but 
at higher resolutions the conclusion is clear: without al- 
gebraic constraint projection the numerical solution does 
not converge to the continuum solution. Even with the 
constraint projection we do not seem to obtain perfect 
second order convergence, but since we are evolving ran- 
dom noise and the trend is clear, we conclude that with 
constraint projection the scheme is converging. Further- 
more in our tests with smooth initial data we do obtain 
perfect second order convergence with the constraint pro- 
jection. We obtain similar results with a Courant factor 
A = 0.1, and with the harmonic slicing and zero shift. 
We also find similar results when discrctizing at fourth 



order. 



Appendix B: Symmetric Hyperbolicity of Z4 with 
the puncture gauge 

Definitions: A quantity E conserved by the principal 
part of the second order system with state vector u, 



E ^ J edx, 

e= {djU,dtu)W^{diU,dtu)\ 



(Bl) 



is called a candidate energy. The Hermitian matrix H is 
called a candidate symmetrizer. A system that admits 
a positive definite candidate symmetrizer is called sym- 
metric hyperbolic. This definition is equivalent to the 
existence of a first order reduction of the system that is 
symmetric hyperbolic according to the standard defini- 
tion for first order systems. 

Z4-C second order form: The principal part of the 
equations of motion of Z4c coupled to the puncture gauge 
can be written 



1 



:apsy 



d? 



'o7y 



a\ lis' 



(B2) 
(B3) 



(B4) 



where we write 9o = [dt — P^di)/a. We thus express the 
principal part of the system in the form 



Sou, 



kl 



(B5) 
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with 

and the principal part matrix given by 



p inn 







kl 



with Iki^^" the appropriate identity and 



kl — 







kl 



P3 






(B6) 



(B7) 



6MS7 7 



ran pk iis_ pk 



T 




where 



Because of the structure of the equations of motion and 
our choice of gauge the block structure of the ansatz can- 
didate is no restriction. The restriction to energy densi- 
ties constructed by using the metric to constract indices 
does restrict the class of symmetrizers, but is the largest 
natural choice. Conservation of the energy is guaranteed 
by Hermicity of the matrix, suppressing non-derivative 
indices, 



for every spatial vector s*, where 

c _ I 
b ,. - I 1 



(B12) 



(B13) 



and the partition of this matrix is compatible with that 
of Uiki into spatial and time derivatives. Imposing energy 
conservation for the puncture gauge implies that 



aPJ mn 
•^\\kl 



(BIO) 



Ansatz candidate and conservation: Starting from the 
ansatz candidate ISOll. 



H 



ij kl nin 



/ TT'i'1 kl nni 


rjij kl 
-"12 











TTi kl m 
-^16 


jjjt ran 
-"12 


-"22 











Tji m 








Tji j k m 
-"33 


Tji k nin 


Tji k 











rrj m kl 
-"34 


Tjkl mn 


Tjkl 











TTjra 
-"35 


-"45 


H55 





, Tjj ran k 
\ -"16 


TT]k 
-"26 











Tjk m 
-"66 



(BU) 



where we define 



Tjij kl ran 
-"11 



IT-i r^M ^,ran , o 7^2 ^M^,k(m^n)l 

+ 2i7fi(7'(''7')-'7'"" -I- 7*(™7")J7'=') 

-"227 7 



H- 



H. 



H 



12 

Trij 
-"22 

ij k m 
33 

i kl ni 



^-^^337 ^ 7 -I- ii.c,a7 ■'7 -t- A337 ^ 7^ 



337 "7 

ini^kl 



337 ' r 



Tji m 
-^26 



h: 



34 

Tji k 



44 

Tjkl 

Tjk m 
-"66 



2i/ie7''(S')™-f i/?g7-^7-, 

^^67™, 

-ff347'^"'7"^'' + -ff|47'S'"", 
-"357 7 

2i7]47'=(™7")' + i/|47'='7"", 
i^i57^-', 

^667""- 



i?l\ 




i^?i 


= 0, 


^?1 


= 0, 


^i\ 


= 0, 


-^12 


= 0, 




= 0, 


-^^16 






_ 1 - MS 


-f^34 


= 0, 


-^44 


= 0, 



-^^34 1 



and the more complicated 



-f^22 — M-L-f^5 



SHl4i2fls - 3) 12H|4(4 + - 4Aii 



24(/25 - 1) 



i2Hl + aHl4), 



1 24 16^25 ^2 , 48(ms-1) 
-^^35 - — 77; 7^^-'^34 -I — — r^-f^44) 



^34 



Ml(4ms - 3/iL)a2 ''''' 

44' 



^1 ^ 2(3/is-4)^2 , 12(^25 -1)„2 



34 "I =2 2 

Ms" 



4^5 - 3/iL 

H 



Hi 



2 _ -v-ps - -;^2 , 12(Ais - 1) „2 



-f^33 - 



2(3/xs - 4) 



fj,sa 



^34 



Note that the the terms antisymmetric in derivative in- 
dices drop out of the conservation equations. 

Positivity: Introducing and expanding with a ten- 
sor basis reveals vanishing elements on the diagonal of 
the candidate symmetrizer; there is no positive defi- 
nite candidate symmetrizer starting from the natural 
ansatz (jBlip . We thus conclude, unfortunately, that the 
use of summation by parts finite differences for the Z4c 
formulation with the puncture gauge will not guarantee 
numerical stability. 
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